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Abstract 

We consider the nonequilibrium work distribution of a quantum oscillator with modulated angular frequency. We 
examine the discrete-to-continuous transition of the distribution as the temperature and the degree of nonadiabaticity 
of the frequency transformation are increased. We further develop a perturbative approach to analyze the effect of 
weak quartic anharmonicities, as well as of a random electric field on a charged oscillator. We find in both cases that 
the degree of nonadiabaticity is enhanced by the perturbation. 

Keywords: nanothermodynamics, nonequilibrium statistics, ion traps 
PACS: 05.30.-d, 05.70.Ln, 05.40.-a 



1. Introduction 

The modern trend of miniaturization leads to the de- 
velopment of smaller and smaller devices, such as nano- 
engines and molecular motors |[T||2][3]. On these very 
short length scales, thermal as well as quantum fluc- 
tuations become important, and usual thermodynamic 
quantities, such as work and heat, acquire a stochastic 
nature. The traditional framework of thermodynamics, 
which neglects fluctuations, thus fails to provide a com- 
plete description of nanosy stems. Extensions of the sec- 
ond law to these small systems have been recently intro- 
duced in the form of the fluctuation theorem ]4] [5) and 
the Jarzynski equality [6|. The Jarzynski work relation, 

<expH?W)> = exp(-/?AF) (1) 

permits determination of an equilibrium free energy dif- 
ference AF from the fluctuations of the nonequilibrium 
work W done on the system during an arbitrary trans- 
formation, quasistatic or not. The classical Jarzynski 
equality ([T} is valid for both isolated and open systems 
Q. It should be noted that the system is initially as- 
sumed to be in a thermal state with inverse temperature 
/?. However, it is not required to be in an equilibrium 
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state at the end of the transformation; for a general pro- 
cess, the system may be arbitrarily far from equilibrium. 
For the case of an open system, the final free energy is 
that of the asymptotic state reached by the system after 
equilibration with the heat bath, while for an isolated 
system, it corresponds to that of the hypothetical equi- 
librium state that would be obtained after coupling to 
the bath. In Eq. ([T} the average of the exponentiated to- 
tal work is taken over an ensemble of many realizations 
of the same process, the nonequilibrium fluctuations of 
the work being characterized by the probability distri- 
bution ViW). The classical Jarzynski equality has been 
verified experimentally using stretched biopolymers (SJ, 
a mechanical torsion pendulum [9|, and a colloidal par- 
ticle in an anharmonic trap IflOl . At the quantum level, 
the Jarzynski relation has been shown to hold for iso- 
lated tlTTirrZl as well as for open systems ifPTlfBl . but an 
experimental investigation is still lacking. This is partly 
due to the fact that quantum work is not an observable 
in the usual sense, as it is not described by a Hermi- 
tian operator in Hilbert space, but by a two-time correla- 
tion function lfT31 . A measurement scheme to determine 
the full quantum work statistics P(W), and hence verify 
the quantum Jarzynski equality, using ultracold trapped 
ions has been proposed in Ref. lfT6ll . We have recently 
explicitly evaluated the continuous envelop of the work 
distribution P(W) for a time-dependent harmonic oscil- 
lator, an analytically solvable model ifTTl . In the gen- 
eral quantum case, however, the work distribution is dis- 
crete, reflecting the quantized nature of the energy spec- 
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trum. In the present article, we provide a detailed dis- 
cussion of the transition from discrete to continuous dis- 
tributions by introducing the cumulative work function. 
We further analyze the effect of a weak anharmonicity 
on P(W) by using time-dependent perturbation theory. 
We finally examine the effect of a weak random electric 
field on the work statistics of a charged harmonic oscil- 
lator, a case of high relevance for the thermodynamic 
study of realistic modulated Paul traps I' TSI . 

2. Time-dependent harmonic oscillator 

We begin by reviewing the solution of the quantum 
parametric oscillator, following the method developed 
by Husimi [19|. We use this opportunity to extend the 
results of Ref. [19| and correct some misprints. The 
Hamiltonian of a quantum mechanical harmonic oscil- 
lator with time-dependent angular frequency to, is 



of the angular frequency a>,. The general form of the 
propagator can be determined from tf/,(x) by noting that 

ifft(x) = J dx U Uo (x\x )tff t0 (x Q ) . (8) 
It is explicitly given by fl9l 



H, 



M 



2M 



2 2 

to.x 



(2) 



The parameterization to, starts at initial value o»o at 
t — to and ends at final value wj at f = t. We 
denote by <f>' n the instantaneous eigenfunctions and by 
E' n = hto, (n + 1 /2) the instantaneous eigenvalues of the 
quadratic Hamiltonian pi. The dynamics of the har- 
monic oscillator is Gaussian for any to,. By introducing 
the Gaussian wave function ansatz, 



M / iM r • 7 7i\ 
U '*= V2« eXP fe[*' X ~ 2XX0 + Y ^\)' 

(9) 

where X, and Y, are solutions of Eq. |7]) satisfying the 
boundary conditions Xq = 0, Xq = 1 an d ^0 = 1> = 0, 
the latter being an expression of the commutation rela- 
tion between position and momentum. 

2.1. Method of generating functions 

The time variation of the angular frequency (|2]i in- 
duces transitions between different energy eigenstates 
of the oscillator. We are thus interested in the transition 
probabilities p T m n from an initial state |n) at to = to 
a final state |m) at t = t. In the following, we use the 
method of generating functions to evaluate pj nn |[T9l . 
We start with the definition, 



Pm,n 



J dx J dx<f>* T m (x)U T , (x\x ) $(x ) 



, (10) 



and denote the complex conjugate of a number z by z*. 
The quadratic generating function of (j)' n {x) is 



if/,(x) 



exp(^\a,x 2 +2b l x + c t \), (3) ^ t r 



the Schrodinger equation for the oscillator can be re- 
duced to a system of three coupled differential equations 
for the time-dependent coefficients a,, b, and c,, 



exp 
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which can be calculated by a Fourier expansion of the 



The nonlinear equation Q is of the Riccati type and is 
therefore solvable. It can be mapped to the equation of 
motion of a classical time-dependent harmonic oscilla- 
tor via a, = MX, IX,, and we obtain 



df- 



X, + to 2 X, = 



(7) 



With the solutions of (|4|-(|7| the Gaussian wave function 
t^,(x) ^ is fully characterized by the time-dependence 



left-hand side of Eq. (Hi. The generating function p T m n 
is then defined as 

P( M ,v) = ^ M "V"/4„. (12) 

m,n 

Combining Eqs. ( [T0| and ([TTJ, we find 

P(u, v) = 

V2 (13) 

Ve*(l - « 2 )(1 - V 2 ) + (1 + M 2 )(l + V 2 ) - AUV ' 

The (m, v)-dependence of the generating function P(u, v) 
remains the same for all possible transformations a>,. 
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Details about the specific parameterization of the an- 
gular frequency only enter through different numerical 
values of the parameter Q* given by, 



Q* = 



1 



[<4 [to] X 2 + X 2 ] + [oj\ Y 2 + Y 2 ]} . (14) 



From a physical point of view, Q* can be regarded as a 
measure of the degree of adiabaticity of the process and 
will be discussed in more detail in the following sub- 
section. Among the properties of the generating func- 
tion P(u, v), it is worth mentioning that the law of total 
probability £ p T m n — 1 is fulfilled and is equivalent to 



p(u, 1) = 



1 



1 — u 



(15) 



For a constant frequency, a> t = wq, we note that the 
solution of Eq. (|7]i is given by 



X, — — sin (u>Q t) , and Y, = cos (coo t) . (16) 
o» 



The latter imply with Eq. ([T4J| that Q* 
thus simplifies to 



P(u,v)\ 



1 



1 — uv 



(17) 



which is equivalent to pj nn - <5 m „, indicating the ab- 
sence of transitions, as expected. The symmetry relation 



of the generating function (13 1, P{— u, — v) = P(u,v), 
further shows that p T m n — if m, n are of different par- 
ity. This is an expression of a selection rule m - n + 2k, 
where k is an integer. In subsection 5.1 we rederive 
this selection rule by means of time-dependent pertur- 
bation theory. We mention in addition that the transi- 
tion probabilities are symmetric, pj nn = p T nnv following 
P(u, v) = P(v, u). Explicit expressions for the transitions 
probabilities p] n „ are given in terms of hypergeometric 
functions in Appendix B 



2.2. Measure of adiabaticity 

The parameter Q* defined in ( 14 1 can be given a sim- 
ple physical meaning |fT9ll . We base our discussion of 
adiabaticity on the equivalent classical harmonic oscil- 
lator |7]i since the generating function P(u, v) ( fT3| > is 
fully determined through the classical solutions X, and 
Y,. For an adiabatic transformation, the action of the 
oscillator, given by the ratio of the energy to the angu- 
lar frequency, is a time-independent constant. For qua- 
sistatic processes we have the two adiabatic invariants, 



X} + co 2 X 2 1 



Yf + uj 2 t Y 2 



and 



CO, 



co 



CO, 



= ojo. (18) 



From the definition (14i of the parameter Q*, we see 
that in this case we simply have Q* — 1 . As mentioned 
earlier this implies P(u, v) = (1 - uv)~ l and pj n n = 5„ h „. 
The latter is an expression of the quantum adiabatic the- 
orem: For infinitely slow transformations no transitions 
between different quantum states occur. For fast trans- 
formations, on the other hand, we can regard Q* as a 
measure of the degree of nonadiabaticity of the process. 
As an illustration, we evaluate mean and variance of the 
energy of the oscillator at time r and express them as 
a function of Q*. For a transition from initial state \n) 
to final state \m), the mean-quantum number of the final 
state (m)„ can be obtained by taking the first derivative 



of the generating function <(T3J> of p T m n , 
Y_u n Yjnp T mn 



dP(u, v) 



dv 

Q*(\ + u) 



(1-H) 



(19) 



2(1 - m) 2 



and expanding the left hand side of ( 19 1 in powers of u: 



1 and Eq. £3} {m)n = £ m ^ = L + I j Q* _ I . 



(20) 



Noting that p\ = exp(-/3E®)/Zo, the mean energy of the 
oscillator at time t then reads 



(H T ) = £ Hcoi {{m) n + i j p° n 
Q* coth^Swo)- 



(21) 



(H T ) 



Q* 



(22) 



Since (m)„ > 0, and hence (H T ) > Hcj\/2, the parameter 
Q* necessarily satisfies Q* > 1 for generic processes. In 
the zero temperature limit, Eq. pT) reduces to 

ha>\ 
~Y 

The above equation corrects a misprint appearing in 
Eq. (5.21) of Ref. Ifl9ll (ojq should be replaced by u)\). 
The mean-square quantum number (m 2 } - (m) 2 at time 
r can be calculated in a similar way by considering 
{mini - !))„. By differentiating Eq. (fT3l twice, we have 



^ M "^m(m-1)^„ 



1 + (u - 6)u + 3Q* 2 (u + I) 2 + 4Q*(u 2 - 1) 
4(m- l) 3 

and a series expansion in powers of u leads to 



(23) 



< m(m - 1) >„ =- [l - 2n (n + 1) - 4 Q* (2n + 1) 
+3 g* 2 (2« 2 + 2« + l)]. 



(24) 
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The mean-square quantum number is then obtained by 



combining Eqs. (20 1 and ( |24| : 

(m 2 ) n - {m) 2 n = ]- (g* 2 - l) (n 2 + n + l) . (25) 



From Eqs. (20 1, (21 1 and (25 I, we can finally write the 



variance of the energy as 



(h 2 )-{H t ) 2 = -^icsch 2 ^ too) 

x (l - 4Q" 2 - 3 cosh ifihojo) + 4g* sinh O8to )) . 

(26) 



The zero-temperature limit, 



(hI)-{h t ) 2 = ^ (e* 2 -i) , 



(27) 



is again the correct version of Eq. (5.22) of Ref. |fT9ll . 
Equation ( p5| indicates that the parameter Q* directly 
controls the magnitude of the variance of the occupa- 
tion number, yn 2 } - {m) 2 . In the adiabatic limit, where 

Q* — 1, we readily get (m)„ = n and (m 2 ^ - {m) 2 = 0. 
We therefore recover that for adiabatic processes the 
system remains in its initial state, \m) - \n). On the 
other hand, for fast nonadiabatic processes, the mean 
(m) n and the dispersion ('« 2 ) - {m) 2 increase with in- 
creasing values of Q*, indicating that the quantum oscil- 
lator ends in a final state \m) which is farther and farther 
away from the initial state |n). The latter corresponds 
to larger and larger final values of the mean energy and 
energy variance, Eqs. ( |2"Tj ) and |26|. 

It is worth mentioning that the above discussion of 
the adiabaticity parameter Q* for the parametric oscil- 
lator is close in spirit to the Einstein criteria for adia- 
batic processes [21 1. Einstein noted that for an adiabatic 
process, the classical action of the oscillator, (H,) /to t , 
should remain constant and the number of quanta should 
therefore remain unchanged. In the present situation, 
we have (H T ) /to T oc Q*, and the action only remains 
constant when Q* — 1 . The latter is precisely the con- 
dition that we derived for an adiabatic transformation. 
For nonadiabatic processes, the parameter Q* > 1 thus 
gives a measure for the increase of the classical action of 
the oscillator. Further discussions of adiabatic measures 
can be found in Ref. [201. 



3. Work probability density function 

In this section, we introduce the probability distri- 
bution P(W) on the nonequilibrium work done on the 



parametric oscillator during a variation of its angular 
frequency. We give the expressions of its continuous 
envelop in the limit of high and low temperatures, for 
adiabatic and nonadiabatic transformations [17|. We 
further compare our nonadiabatic results with those re- 
cently derived by van Zon and coworkers in Ref. [23 1. 

In quantum mechanics, the probability density func- 
tion of work is obtained by considering the difference 
between final and initial energy eigenstates, EJ n and E®, 
averaged over all possible final and initial states. The 
probability distribution of the total work done during 
time t can thus be written as [ 15 1, 

PW = £ 6 (W - {El - E°„)) p T „ hn pi (28) 

m,n 

where p° n = exp(- y 6£ , J ] ) )/Zo is the initial (thermal) occu- 
pation probability. Equation ( |2"8"j ) show that work is a 
random quantity because of the presence of both ther- 
mal and quantum uncertainties, encoded respectively in 
p\ and p T mn . The characteristic function of the work, 
defined as the Fourier transform of the probability dis- 
tribution d28l, 



I 



dWexp(ip.W)P(W), 



(29) 



can be written in closed form in terms of the energies, 
fio = too, fii = toi (Ae - S\ - so). Using expression 
([T3]l of the generating function, we have 



G(p) = J - ex P ( -M exp [in y j (30) 

where the denominator is given by 

n (Q*) =Q* (1 - exp (2^£i)) (1 - exp (-2(ifi +j8)e )) 
+ (1 + exp (2ifisi)) (1 + exp (-2(i/j +y3)e )) 
-4 exp {inei ) exp (-(ip + /3)s ). 

(31) 

The above equations for G(p) are exact and fully char- 
acterizes the work distribution of the time-dependent 
quantum harmonic oscillator Q for arbitrary parame- 
terizations of the angular frequency a> t . As mentioned 
previously, different realizations of at, merely lead to 
different values of the parameter Q*. The direct an- 
alytic evaluation of the nonequilibrium work distribu- 



tion V(W) by Fourier inverting Eq. (30 1 is not feasible 



in the general case due to the nonanalytic denomina- 
tor of G(p). We here provide analytical approximations 
in various limits of interest. The detailed derivation of 
the following probability distributions can be found in 

Ref. na. 
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Adiabatic transformations. As discussed in subsection 
2.2 adiabatic transformations are characterized by Q* = 
1. By approximating Eq. ( |30] ) in the limit of zero- 
temperature, h/3 » 1, we get 



P(W) = 6 \W- 



Ae 



(32) 



Equation ( 32 1 expresses the deterministic nature of adia- 
batic processes at zero-temperature: the oscillator starts 
and ends in its ground state. Work is hence simply given 
by the difference of final and initial ground state ener- 
gies. In the opposite classical limit, h/3 <c 1, we find 

PiW) = &l ex J^ W ) &m (33) 
Aco \ Aco I 



Equation ( 33 1 is identical to the classical work probabil- 



ity distribution derived by Jarzynski E2l . Note that in 
this case only positive work fluctuations occur, whose 
magnitude is controlled by the finite temperature. 

Nonadiabatic transformations. We next consider nona- 
diabatic transformation which correspond to Q* > 1. In 
the zero-temperature limit, hfi » 1, the work distribu- 
tion can be approximated in the limit of small E\ by, 



exp(- 



W-As/2 \ 



>MG*-l)£i(W-A£/2) 



(34) 



The zero-temperature nonadiabatic work distribution 
p4| ) is valid when W > As/2. This condition follows 
from the existence of the minimal ground state energy 



of the oscillator. An expansion of Eq. ( 30 1 in the classi- 
cal limit, h/3 <K 1, leads to 



/3 2 co 2 /n 2 



y 2Q*cl>qlo\ - cJq — co 2 
Q*co coi - col 



x exp 



xK 



2Q*co coi - col ~ w i 

ptoQcox Ve* 2 - 1 

\2Q*LOocoi - Wg - co\\ 



/3W 



\W\ 



(35) 



where T(x) denotes the Euler Gamma function and 
K v (x) is the Macdonald function, that is the modified 
Bessel function of the third kind. Equation ( (35| ) is well- 
defined provided the term under the square-root is posi- 
tive. This implies that the parameter Q* is larger than 
the value (2* s that corresponds to a sudden (instanta- 
neous) switch of the frequency from coq to co\ , 



9 9 
cot + CO, 

q* > ^ — - = e: s • 



2co()CO 



(36) 



<KW) 

l.Or 




w 



25 



Figure 1: Cumulative work distribution <j>(W) (38) for the exact quan- 
tum expression (red) and the approximated classical form J35J 
(blue, dashed) for the parameters tog = 1, b)\ = 1.5, h = 1,/? = 0.5 
and Q* = 3. The discrete structure of the quantum distribution is eas- 
ily recognized as compared to the continuous classical distribution. 



An example for a frequency parameterization resulting 
in an unbounded, divergent value of Q* can be found in 
Ref. J24). In the regime where Q* < g* s , the work dis- 
tribution has been derived by van Zon and coworkers by 
explicitly evaluating the complex integral of the inverse 
Fourier transform ll23l. It reads in our notation 



<P(W) 



cof, — 2Q*coqco\ + co 2 . 



x exp 



x/ 



Q*coqco x - co\ 

^■Q*CO()COi — CO 2 — Cl 

\2Q*co coi - co 2 - co 2 



-pw 



(37) 



®(W), 



where Iq(x) denotes the modified Bessel function of the 
first kind. It is worth mentioning that, in contrast to 
the adiabatic ( |33| l and nonadiabatic ( |3"7| ) results, nega- 
tive work values can occur in Eq. ( [35] l for large values 
of the parameter Q* 11231 . 



4. Transition between discrete and continuous 
nonequilibrium work distributions 

We now turn to the analysis of the transition from 
the discrete work probability density (|28|i to the con- 



tinuous analytical approximations Eqs. (33 I, (34i, (35i 



and ( 37 1, and provide a criterion for the observation of 

The definition 



the discreteness of V(W). The definition ( |28| > makes 
clear that the nonequilibrium work distribution consists 
of a discrete sum of delta peaks when the energy spec- 
trum of the system is quantized. This is in stark con- 
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trast to classical work distributions which are continu- 
ous. The discreteness of P(W) can therefore be regarded 
as the hallmark of the quantum nature of work. In order 
to quantitatively investigate the discrete-to-continuous 
transition, we introduce the cumulative probability dis- 
tribution <p(W) defined as, 



<f>(W) 



dW P(W) 



(38) 



In the above equation, the constant x nas to be chosen 
according to the range validity of the work distribution, 
e.g. in the zero-temperature limits = As/2. The cu- 
mulative distribution ( |3~8"| > is a staircase function in the 
case of a discrete distribution P(W) and turns over to a 
smooth function in the continuous limit. Figure[TJshows 
the exact quantum function 4>(W) for a given set of pa- 
rameters, together with the continuous classical approx- 



imation corresponding to Eq. ( 35 1. The numerical value 
of the parameter Q* has been chosen to describe realis- 
tic experiments with ion traps |26l . In order to obtain 
a criterion for the discreteness of the work distribution 
ViW), we consider the mean energy of the harmonic os- 
cillator as given by Eq. (21 1. The energy spectrum of the 



harmonic oscillator is usually considered to be quasi- 
continuous when the mean quantum number is much 
larger than the level separation |25l . or in other words, 
when the mean energy is much larger than the energy 
quantum, (H T ) » Hu>\. An inspection of Eq. (21 1 re- 
veals that this limit can be achieved in two ways: 

High temperatures. In the limit of high temperatures, 
hp «: 1, we have coth(/3/2 Hojq) » 1. For a fixed value 



of Q*, we therefore expect the nonequilibrium work dis- 
tribution to become continuous when the temperature is 
increased. Figure [2] depicts the exact cumulative distri- 
bution 4>{W) and the approximate result corresponding 
to Eq. ( 35 1 for the same parameters as in Fig. [TJexcept 
the temperature which has been increased by a factor 
ten; the two curves are indistinguishable. 

Large Q* values. Another way to reach a large mean 
energy pi) is to increase the value of Q*, while keep- 
ing the temperature constant. The work distribution can 
thus become continuous when the degree of nonadia- 
baticity of the process is increased. Figure [3] shows the 
exact and approximate cumulative distribution 4>(W) for 
a value of Q* twenty times larger than in Fig.[T] all other 
parameters being the same: the two curves are again not 
distinguishable. It is interesting to notice that the tran- 
sition to a continuous work distribution occurs faster 
when augmenting the temperature than the degree of 
nonadiabaticity Q* (a factor two for the parameters of 
Fig.|TJ, the reason being that the mean energy pT) does 
not depend linearly on temperature in contrast to Q* (the 
temperature dependence becomes linear only for high 
temperatures). 

Noninteger ratios si/bq. One may wonder whether it is 
possible to have continuous work distributions, for any 
values of Q* and by for instance considering non- 
rational quotients s^/sq. In this limit, the gaps between 
d>-peaks could be filled and V(W) would appear contin- 
uous. This is not the case, however. We note indeed 
that for any ratio ei/eo, we do not obtain real contin- 
uous distributions with respect to the real set. The 5- 



functions in Eq. (28 i only contribute to the total sum if 
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w 



Figure 4: Cumulative work distribution (38) for an initially pure state 
|«o)> illustrating the sum rule, m = n±2k, which limits the possible 
sizes of the work steps (;iq = 3, £q = 1, £i = 1.5 and Q* = 3). 



their argument is zero. The only permitted work values 
are therefore of the form, 



W = £\\m + 



1 



s \n + 



1 



(39) 



Now since m and n are integers, only countably infinite 
values of W can occur, implying that the work distribu- 
tion is discrete in any case. The latter follows from the 
fact that rational numbers are densely distributed within 
the real set. 

Effect of the selection rule. Finally, we examine the ef- 
fect of the selection rule, m — n ± 2k, noted in sec- 
tion[2]on the work distribution P(W). The existence of 
this selection rule limits the possible values that the step 
sizes of the cumulative distribution <p(W) can take. For 
simplicity, we consider a harmonic oscillator which is 
initially prepared in a given energy eigenstate |no). A 
direct consequence of Eq. < [39] > and of the selection rule 
is then that the allowed work values are, 



W = h (0l>\ - Wq) «o + 



1 



hu> x 2k, (40) 



where k are integers. For two neighboring work values, 
the step size is constant and reads, 



AW = 2h(jj x 



(41) 



In Fig.|4]we illustrate the effect of the selection rule by 
plotting the corresponding cumulative distribution with 
the same parameters as in Fig.[T] 

5. Perturbation theory 

We develop in the present section a perturbative 
approach to determine the quantum work distribution 



P(W). We treat in detail the case of a small anharmonic 
correction to the potential, as well as the effect of a small 
external fluctuating electric field on a charged harmonic 
oscillator. Both situations are motivated by the exper- 
imental study of the quantum work statistics in linear 
Paul traps [16|. We begin by giving a simple perturba- 
tive derivation of the sum rule of section [2] which does 
not require the full solution of the Schrodinger equation. 



5.7. Selection rule 

We start by rewriting Hamiltonian |2} as 

H t = Ho + £l t , 



(42) 



where we have introduced the unperturbed Hamiltonian 

(43) 



tt P 2 ^ M 2 2 
tin — U) n X 

2M 2 



{<4 - ^) 



(44) 



and a "perturbation" term 
M 

"T 

The latter can be considered small for small frequency 
changes. In first order time-dependent perturbation the- 
ory, the transition probabilities between initial state |n) 
and final state |m) are given by ll27l . 

2 



5 m ,n + j dt exp(/w„,,„r) Qj„ „ 



(45) 



where Ha>„ 



EZ - £2 denotes the difference 



of the unperturbed energy eigenvalues and QJ nn = 
{m\Q. t \n) are the corresponding interaction matrix el- 
ements. By expressing the position operator, x = 
y]hl2Mu)Q (a? + aj, in terms of the usual ladder oper- 
ators, cft\n) — Vn + 1 |n +1) and a \n) = yfn\n - 1), the 
interaction matrix elements can be written explicitly as 



Q'm n = ~- A — ( w o - [ Vn+ 1 in + 2 (5 m> „ +2 
+ (2n+ 1) S m j, + Vn - 1 Vn5„,,„_ 2 ], 



(46) 



with the Kronecker-delta 6 m n . Equation ( |46| > shows that 
only transitions that satisfy m - n + 2 are possible, 
which is precisely the selection rule noted earlier. It 
should be emphasized that this selection rule is at vari- 
ance with usual textbook examples which contain the 
selection rule m — n + 1 . The latter applies to a quan- 
tum oscillator driven by a small perturbation linear in 
the position, whereas we here deal with a perturbation 
( |4"4"1 > which is quadratic in x. The full expression of the 
transitions probabilities (|4"5j) that follow from Eq. (|4"6]> 



is given in Appendix C 
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1.00 1.01 1.02 1.03 1.04 1.05cm 



Figure 5: Variance and mean (inset) work for an oscillator with weak 
anharmonic corrections (red) j47| compared with those of the unper- 
turbed oscillator (blue, dashed) |2J (a> = 0.5, r = 1,jS = 0.5, M = 1, 
h = 1 and a a = 0.025). 



5.2. Anharmonic corrections 

A method to experimentally measure the quantum 
work distribution in modulated ion trap systems has 
been proposed in Ref. ifTBI . In these systems, the confin- 
ing potential is harmonic to a very good accuracy ||26l . 
One attractive feature of linear Paul traps is however the 
possibility to modify the shape of the potential with the 
help of external gate voltages. We here investigate the 
influence of a small quartic anharmonicity on the work 
distribution P(W). As before, we write the total Hamil- 
tonian as 

H, = H + Q, + A, , (47) 



where the first anharmonic correction is given by, 



A, 



(48) 



The total transition probabilities can then be written as 

T 2 



fofdt exp (iecv, f) (n;, v , + A' m n ) 



(49) 

where A' mn are the anharmonic interaction matrix el- 
ements. The analytic transition probabilities p T mn are 
again given in |Appendix C| For the sake of clarity, we 
will continue with a numerical discussion of the results. 
For the numerical analysis, we choose the parameteri- 
zation of ajf to be linear in time, 



a>, - oj,, + 



(50) 



Since the anharmonic corrections are given by the geo- 
metric set-up and hence directly scale with the angular 
frequency of the harmonic oscillator, we assume that 



a, = o- a [u\ - «q) - . 



The parameter cr a controls the strength of the perturba- 
tion. In Fig. [5] we have plotted the mean work (W) 
and the variance, var(W) = (w 2 ) - (W) 2 , of the work 
distribution for the anharmonically perturbed harmonic 



oscillator ( 47 1, together with the exact result for the un- 



perturbed oscillator pt. We observe that both quantities 
are enhanced by the anharmonicity. From the analyti- 



cal expressions of the transition probabilities (C.3 I, we 
see that additional transitions, m — n + 4, now become 
possible because of the quartic correction. These ad- 
ditional transitions lead to a larger mean and variance 
of the work. Based on our discussion in section |2~2"1 we 
can therefore conclude that the anharmonic perturbation 
increases the degree of nonadiabaticity of the frequency 
change. Numerical comparison further shows that the 
effect of A, can be neglected up to a strength of roughly 
one percent, cr a < 0.01, of the harmonic amplitude Q,. 
For a standard trap configuration with trap frequencies 
of the order kHz-MHz, the harmonic assumptions is ful- 
filled up to energies of the order of eV, see Ref. [26 1, 
and the effect of anharmonic corrections are negligible 
for these energies. 

5.3. Random electric field corrections 

Linear Paul traps are almost perfectly isolated from 
their surroundings. They however suffer from the pres- 
ence of random electric fields that are generated in the 
trap electrodes [28 1. These weak fluctuating fields are 
the source of motional heating of the charged ions con- 
fined in the harmonic trap. The Hamiltonian of the 
quantum oscillator in the presence of the field is 

H, = H Q + n, + A, (52) 

where the small perturbation A, is linear in position, 

A, = A t x . (53) 

The function A, = qE, is proportional to the random 
electric field E, (q is the charge of the ion) and is taken 
to be Gaussian distributed with 



(A t ) = and (A, A s ) = K tl 



(54) 



(51) 



The heating rate of the trap is related to the spectral den- 
sity of the noise A t l29l 



<») — T~rrz f ds exp(iu),s) (A,A t+s ) . (55) 
A-Mhu), J 



We first calculate the transition probabilities p T m n for a 
fixed value of A, and then average over A, using Eq. (j54j. 
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In complete analogy to Eq. (49 1, we obtain 

r 

Pm,n = 6 >n,n + f dt eX P « f ) ( Q m,« + A m,J 


with the interaction matrix elements A r m n given by, 
A, 



(56) 



A' 



^2M 



(V«+ 1 1 + V«ev»-i) 



(57) 

The explicit expression of the transition probabilities 



can be found in Appendix C.3 After averaging over all 
possible A,, the transition probabilities can be divided 
into two distinct contributions coming from the para- 
metric variation of the frequency (£2, in Eq. ((52])) and 
the noise term (A, in Eq. ((52])), 

(?»,,»),,, = PlnjM) + Pm,n « A > **» • ( 58 > 

Similarly, we can separate the mean final energy into a 
deterministic and a stochastic part, 



(H T ) = J] E] n (p T m , n (a> t ) + pl hn ((A, A s )j) p° n 

(e* + <2i,)coth(| ruo o y 



Here the parameter Q*. is defined as 



(59) 



fuoi/2 coth(/?/ '2 hcoo) 



(60) 



with 



(H) M = J] ^ ( m + Pm,n «^ A*» />« • (61) 



Equation ( |49[ ) shows that the effect of the random elec- 
tric field is to renormalize the adiabaticity parameter 
Q* — > Q* + Q* A . Both the mean and the variance of 
the work distribution are increased as depicted in Fig. [6] 
The fluctuating field thus enhances the degree of nona- 
diabaticity. This effect can be understood by noting that 
the perturbation A, generates additional transitions be- 
tween states (the latter obey m — n ± 1). We observe 
that the variance is more sensitive to the perturbation 
than the mean, since it depends quadratically on Q* and 
not linearly. For the numerical analysis we have chosen 
a white noise of the form, 



K t , s = cr A \cd>l - cjfj 5(t - s) , 



(62) 



where the relative noise strength is given by cr A . As for 
the case of the anharmonic perturbation, we note that 
one can neglect the influence of the electric noise up to 
a relative strength of roughly one percent, cr A < 0.01. 



var(W) 
0.015 

0.010 

0.005 
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Figure 6: Variance and mean (inset) work for a charged oscillator with 
weak electric noise (red) j47) compared with those of the unperturbed 
oscillator (blue, dashed) {2} (w = 0.5, t = l,/3 = 0.5, M = 1, ft = 1 
and cr„ = 0.025). 



6. Conclusion 



We have considered the statistics of the nonequilib- 
rium work of a quantum oscillator when its angular fre- 
quency is varied in time. Due to the quantized nature of 
the energy spectrum, the work probability distribution is 
discrete. We have analyzed the discrete-to-continuous 
transition of the work distribution in various limits by 
introducing the cumulative function. We have shown 
that the cumulative work distribution becomes smooth 
in the limit of high temperatures and of large values of 
the adiabaticity parameter Q*; in both regimes, the mean 
energy of the oscillator is much larger than the energy 
separation, and the spectrum can be considered quasi- 
continuous. We have moreover developed a perturbative 
approach to investigate the effects of small quartic an- 
harmonicities on the work distribution. We have found 
that the latter increase both the mean and the variance of 
the final energy of the oscillator, indicating an augmen- 
tation of the nonadiabaticity of the frequency change. In 
a similar way, we have studied the influence of a weak 
electric noise on a charged harmonic oscillator and ob- 
tained an analogous enhancement of the degree of nona- 
diabaticity. Our results permit an accurate description 
of measured quantum work distributions in modulated 
Paul trap under realistic experimental conditions. 
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Appendix A. Analytical expression for Q* 



Appendix B. Exact transition probabilities 



Closed expressions for the adiabaticity parameter Q* 
can be found whenever the classical equation (|7]i 



dt 



2 x, + co;x t = o 



(A.l) 



can be solved analytically. Equation (jATTji is of the gen- 
eral form of a Hill equation which can be solved under 
various conditions [ 30 , 31]. Equation (A.l i reduces to 



common differential equations in the case of specific pa- 
rameterizations u,. Thus for the case of linear parame- 
terizations ( |50"| ), the solutions are given in terms of the 
Airy-functions [ 17 1. On the other hand, for a sinusoidal 
parameterization 



2 2/2 

or, = o), — o>i — to, 



(o>\ - tug) cos (|-) (A.2) 



Eq. ( A. 1 1 takes the form of the Mathieu equation. The 
solutions X, and Y, can then be written as, 



X, =— |S [rfi,rf 2 ,0]C 

7T 



d\,d.2, 



nt 
At 



-C [di,d 2 ,0]S 



d\,dz, 



nt 
A~t 



(A.3) 



X C [£/i,£/ 2 ,0] S [£/!,£/ 2 ,0] 



-1 



and 



-C [d l7 d 2 ,Q]^ [^i,^2,0] 



d\,d 2 , ~r 
At 



Y, =|S' [d u d 2 ,0]C 
-C [d u d 2 , 0]S 
x|C [Ji,d 2 ,0]S' [d u d 2 ,0] 



(A.4) 



di,d 2 , — 
At 



-C [d l ,d 2 ,0]S[d l ,d 2 ,0] 



Here the functions C and S denote the corresponding 
Mathieu functions [32]]. The parameters d\ and d 2 are 
given by, 

16 o T 

di = -z oj 2 ,t 2 (A.5) 



and 



I 2 2\ 2 



(A.6) 



Further analysis of the parameter Q* in the context of 
vacuum squeezing can be found in Ref. |24|. 



We here collect the analytical expressions of the tran- 
sition probabilities p T mn |fT9ll . Despite its apparent sim- 



plicity, the generation function P(u,v) (13i cannot be 



expanded in powers of u and v in an exact series. We 
thus make use of the p] n n as defined by the matrix ele- 
ments of the propagator U Tfi (x\x ), p] n n = \U r m J^ |lo| >- 
This matrix elements are given by 

Ul,n = j dx j dxf T m {x)U % o{x\xo)<p n (xo). (B.l) 

We use again the method of generating functions. We 
use the linear generating function of tp'^x) G71 . 



n=\ x ' 



JMoj, 



exp 



Ma), 



x z + 2 



Mix), 



(B.2) 



■zx-z 



2h V h 

to evaluate the generating function of the propagator, 

^)=zud MV,,t/ -- (B3) 

By introducing the complex parameters, 



(jj\a>aX T - u>()iX T + <jj\iY T + Y T 
Iff = luy^x (Q*- 1) 

cr| 2 = 2w Wl(G* + l) 



(B.4) 
(B.5) 



we can write 



■ya>QLOi I - Ai -dajoujiuv + £*v 
U (u, v) = - exp 



cr 



(B.6) 

The matrix elements U T m n can then be obtained by a se- 



ries expansion of (B.6 1 in powers of u and v [ 19 



<v» -v^^owi y 2 „ +m _ 1 . (r „ + , 



■/i+m+ 1 



mill (m,n) 



4^ l\ \(n 



[-2iV2/(e*-l)]' 



(B.7) 



![(«-0/2]![(m-/)/2]! 



According to the selection rule m - n ± 2k, I runs over 
even numbers only, if m, n are even, and over odd num- 
bers only, if m, n are odd. The explicit expression for 
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the matrix elements UL „ then reads for even elements 



^2/z,2y 



22v+2/.(-l ; - 



'0^+2^+1 TQi+V)T(y+l) 
1 2 



(B.8) 



x 2 Fi \-fi, -v; 



2' 



and for odd elements 



U 



2/i+l,2v+l 



|8i(2v+ l)!(2/i + 1)! 
(Q* - l)2 2v+2 f +1 



2v+l r* 2/1+1 



(B.9) 



o-iv+^+i rou + i)r(v + i) 

3 2 



x 2j Fj -yu, -v; 



2 1-Q* 



We have here introduced the hypergeometric function 
iF\ 11321 in order to simplify the sums and write the 
matrix elements U T m n in closed form. T(x) denotes the 
Euler Gamma function. Combining everything, we get 
the explicit expressions for the transition probabilities 
which reads for even transitions 



2 l/2 



e*-i 



jJ+V 



Pi^y ( q* + ^1/2 \Q* + \ 

r(i/2 + /i)r(i/2 + v) 

nT(l + jU)T(l + v) 

/ 1 2 

2F l -v; -; — 



(B.10) 



and for odd transitions 

2 7/2 

P2u+l,2v+l 



X 



6*-i 

(2* + 1)3/2 + 1 

r(3/2+yu)r(3/2 + v) 

7rr(i +yu)m + V) 



CB.ll) 



2 Fi -v; -; 



2' 1-Q? 



Appendix C. Perturbational transition probabilities 

In this appendix, we provide the approximate first- 
order transition probabilities calculated with the help of 
time-dependent perturbation theory. 



Appendix C.l. Isolated harmonic oscillator 



Following Eqs. (45 i and (46 1, the transition probabil- 
ities for an isolated harmonic oscillator are 



Pm,n 



T 

J dt exp (;w,„,„r) (w„ - w?) 



Vn + 1 Vn + 2<5 m ,„ +2 + (2n + 1) 5 m , n (C1) 



+ -\fn Vn - 1 (5 m> „_2 

The selection rule m = n ± 2 finds its origin in the pres- 
ence of the Rronecker deltas. 

Appendix C.2. Anharmonic corrections 

The interaction matrix element A' m n (48 1 are given by 

Vn + 1 Vn + 2 Vn + 3 Vn +4 £,„ „ +4 



+ (4n + 6) Vn + 1 Vn + 2 (5 m> „ +2 

+ (6n 2 + 6n + 3) <5 mj „ 

+ (4n-2) V«Vn-l^ m ,„-2 

+ Vn Vn - 1 Vn - 2 Vn - 3 <5,„.„_ 4 



(C2) 



The transition probabilities (49 1 for a harmonic oscilla- 
tor with quartic corrections then read 



P m,n 



r 



4tL» 



(wy-w 2 ) Vn + 1 Vn + 2 



+ (2n + 1 )£„,,„ + Vn Vn - 1 (5,„,n-2 
*2 



+a< 



Vn + 1 Vn + 2 Vn + 3 Vn + 4 <5 mj „ +4 

+ (4n + 6) Vn + 1 Vn + 2<5 m „ +2 

+ (6n 2 + 6n + 3^ 5 m> „ 

+ (4n - 2) V" Vn - 1 £„,,«-2 

+ Vn Vn - 1 Vn - 2 Vn - 3 5 m „_4 



(C.3) 

where additional transitions m = n ±4 become possible. 
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Appendix C.3. Random electric field corrections 



In the presence of an external electric field (53 1, the 
transition probabilities (|56]l become 



Pm,n = 



T 

S m ,n + J dt exp (ico m>n t) 



4co 



(oJq - ) Vn + 1 Vn + 2 £ B!i „ + 2 



+ (2n + 1) + Vn Vn - 1 £ m ,„_2 



\ 2Mwo 



+ 15, 



m,«+l 



V»<W-l) 



(C4) 



The Kronecker deltas now also allow next-neighbor 
transitions m = n±\. Averaging over the noise A, using 



Eq. ( 54 1, we finally obtain, 



T 



xj - (w„ - c^) Vn + 1 Vn + 2 £„,,«+ 



+ (2n + 1) + V« Vn - 1 (5 m> „_ 2 
[(n + 1) £„ v , + i +nd m>n -i] 



2M ho) Q 



l i 

^ dt ^ ds exp (z*a> mj „ (f - s)) {A } A s ) 



(C5) 
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